Infrared emission from interstellar dust cloud with two 
embedded sources: IRAS 19181+1349 
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Abstract 



Mid and far infrared maps of many Galactic star forming regions show multiple 
peaks in close proximity , implying more than one embedded energy sources. With 
the aim of understanding such interstellar clouds better, the present study models 
the case of two embedded sources. A radiative transfer scheme has been developed 
to deal with an uniform density dust cloud in a cylindrical geometry, which includes 
isotropic scattering in addition to the emission and absorption processes. This scheme 
has been applied to the Galactic star forming region associated with IRAS 19181+1349, 
which shows observational evidence for two embedded energy sources. Two indepen- 
dent modelling approaches have been adopted, viz., to fit the observed spectral energy 
distribution (SED) best; or to fit the various radial profiles best, as a function of 
wavelength. Both the models imply remarkably similar physical parameters. 
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1 Introduction 



Galactic star forming regions mostly comprise of Young Stellar Objects (YSOs) / protostars 
still buried inside / in the vicinity of the parent interstellar cloud from which they are 
formed. Hence the study of YSOs leads to understanding of the interstellar medium in the 
close neighborhood of the starbirth. Early evolution of star forming regions is obviously more 
important from the point of understanding the star formation process itself. Typically the 
photons from the embedded energy source, protostar / ZAMS star, get reprocessed by the 
dust component of the interstellar medium in the immediate neighborhood. The dust grains 
absorb / scatter the incident radiation, depending on their dielectric properties. The dust 
grains acquire an equilibrium temperature based on the local radiation field which depends 
on the distance from the energy source and secondary radiation from other grains. The 
reradiation from the grains in the outer regions, is what is observable. Hence, the emerging 
observable spectrum can in principle be connected to the spectrum of the embedded protostar 
/ ZAMS star through detailed radiation transfer provided some details about the geometry 
are known. 

In order to study the earliest stage of star formation, a large amount of observational ef- 
fort is directed towards far infrared / sub-mm observations of prospective young star forming 
regions. Many Galactic star forming regions have been mapped with near diffraction limited 
angular resolutions, using Kuiper Airborne Observatory, Infrared Space Observatory (ISO), 
James Clerk Maxwell Telescope etc. in recent times. In many cases, the maps of continuum 
emission resolve several closeby individual intensity peaks as evident from the morphology 
of their isophots. 

The present study is a step in the direction of extracting maximum possible information 
about the geometrical and physical details of the source by comparing radiative transfer 
models with observations in cases where two nearby sources are resolved. Here, the "nearby" 
implies the interference in energetics of each of the two resolved source by the other one. 
The heating of dust in a cloud with multiple sources has been studied by Rouan (1979) 
under certain simplification. However, the problem of more than one embedded source has 
rarely been addressed quantitatively. The existing observational data showing clear evidence 
of resolved multiple embedded sources justify the need to explore geometries dealing with 
more than one embedded source. Ghosh & Tandon (1985) attempted to study the case 
of two embedded sources with many simplifying assumptions. They neglected some basic 
phenomena like scattering which limited its applicability to A > 50 /im only. The present 
work is an extension of this earlier attempt by including the effect of isotropic scattering. 

In section 2 the problem has been formulated and the radiation transfer scheme is 
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described. In section 3, we model observations of IRAS 19181+1349 which shows evidence 
of having two embedded sources. The results of our modelling are then discussed. 

2 Model Formulation 

The primary aim of the present model is to reproduce the infrared emission from a star 
forming cloud with two embedded sources, keeping the computational complexities at a 
minimum level. 

The interstellar cloud is assumed to be of cylindrical shape. As a starting point, an 
uniform density of the cloud has been assumed. The line joining the two embedded ZAMS 
stellar / protostellar energy sources defines the symmetry axis of the problem. Around each 
of the sources there will be a dust free cavity (Fig. 1). The existence of such a cavity 
is widely accepted due to evaporation of the dust grains in the intense radiation field. In 
addition, radiation pressure effects on the dust grains may also play a role in deciding the 
cavity size. The radiation transfer is carried out through the dust component alone. Dust 
grains with a continuous size distribution have been considered, and their composition is 
a parameter of modelling. Three types of grains, viz., Graphite, Astronomical Silicate and 
Silicon Carbide have been invoked since their existence is generally accepted. All properties 
of the dust grains, viz., absorption and scattering coefficients as a function of wavelength, 
for various sizes and all three types of grains, have been taken from Laor & Draine (1993). 
The size distribution of all the three types of grains has been assumed to be power law 
(n(a)da ~ a~ 3 5 da) as per Mathis, Rumpl & Nordsieck (1977). The wavelength grid used 89 
points covering from the Lyman continuum limit to the millimeter wavelengths. 

The geometrical parameters relevant to the model (Fig. 1) include the radius of the 
cylinder (R C yi), radii of dust free cavities near the two sources ( R c i, R C 2), and distance 
between two exciting sources (D). Other physical parameters are the composition of the 
dust (relative abundances of three types of dust components) and the dust density expressed 
in optical depth at 100 /zm (r W o). The optical depth at any other wavelength is uniquely 
connected to Tioo via the dust properties and composition assumed in the perticular model. 
The cylindrical interstellar cloud is divided into n z identical discs. Each of these discs are 
further divided into n r annular rings. The n z is chosen such that each disc is optically thin 
even at UV, along the z-axis. Along the radial direction (total number of grid points being 
n r ), a two stage grid has been employed which is initially nearly logarithmically spaced (near 
the symmetry axis) and linearly spaced in the outer regions of the cylinder. This scheme 
of radial grid has been arrived at by keeping the optical depth related inaccuracies under 
check, for the entire wavelength region considered for the radiation transfer. Both the near 
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logarithmic and the linear grids are matched by ensuring that radial cell size, 5r(n r ), is a 
smooth function of n r . For modelling attempts of IRAS 19181+1349 a grid of 600 points in 
axial direction and of 25 points in radial direction were employed. 

The relevant calculations are represented in equations [l|-|7|. For clarity the frequency 
suffix has been dropped from all the terms, though calculations are performed for each 
frequency grid point. The quantities in angled brackets are averages over the dust size 
distribution. The code is simplified and optimized in many ways in view of the memory 
requirements and speed. Initially, factors totally dependent on the geometry of the problem 
and which do not change in each iteration are calculated. These include the optical depth 
terms and the geometric integrals involved in computations of radiation received by an unit 
volume element of the ring i from the ring j (equation [1]). The geometric symmetry is 
such that such terms only depend on the axial separation between the two rings and their 
respective radii. The total flux absorbed and scattered by unit volume in each ring due to 
radiation from other rings is then calculated (equation |2[). Also, the radiation received by an 
unit volume element of the ring i from embedded sources depends only on the geometry and 
the optical depth per unit length of the particular model and hence is fixed over iterations 
(equation 0) . The equilibrium temperature of the dust grains (on the median circle) for any 
particular annular ring is calculated using an iterative scheme by equating the power radiated 
by the dust (equation |) to the power absorbed (equation 0). The latter is contributed by the 
embedded exciting sources (attenuated by the line of sight dust) as well as secondary emission 
and scattering from dust grains in all other annular rings (equation |J). This simplifies 
the calculations leading to a set of coupled equations with only two parameters per ring, 
temperature ( Ti ) and F l changing from iteration to iteration thus greatly reducing the 
memory requirements . 

K' J = E ~ T2 [< ™*Q<*» > n « B ^) + A TA )ArjAzA9rj (1) 

V d h,e An < ™ 2 Q ext > 

n r xn z 

Fjl= £ < 7ra 2 Q ext > n d F? (2) 
F i = j2< na 2 Q ext > ^El e -<^Q^>n d {n, s -r c , a ) ^ 

s=l T i,s 

F i = F; + Fl (4) 
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Fl =< ira 2 Q abs > AnB{Ti) 



(5) 



PEmitted = J F\ dv ( 6 ) 
P Absorbed = [ ~ Sf"* ~ (F* + F % )dv (7) 

J < na z Q ext > 

where 

d 2 = ( Zi - Zj ) 2 + (n) 2 + (r,) 2 - 2nrj cos(#) 

9 is azimuthal angular difference between two volume units under consideration. 

r ls = (zi-z s ) 2 + rf 

r 2 s = dust cavity radius for source s. 

Qext = Qabs Q scat 

n d = number density of dust grains. 

Ti = Temperature of dust in ring i. 

F s = surface flux spectrum for source s. 

F l r = Flux absorbed and scattered by unit volume of ring i, due to other rings. 
F l s = Flux absorbed and scattered by unit volume of ring i, due to sources. 
F % = Total flux absorbed and scattered by unit volume of ring i 
-B(T) = Planck function. 

The dust temperature for each annular ring, leading to the temperature distribution 
throughout the cloud, is determined iteratively. Initially (the very first iteration), only the 
two embedded sources power the heating of the grains. From the second iteration onwards 
the effects of secondary heating and scattering are taken into account. In this iterative pro- 
cedure the temperature of each annular ring is gradually updated in each iteration satisfying 
the condition PEmitted = PAbsorbed- The iterations are continued till the fractional changes 
in absorbed power for each annular ring, between successive iterations, reduces below the 
convergence criteria. 

The emergent intensity distribution, as seen by a distant observer is predicted by inte- 
grating the emitted and scattered radiation along relevant lines of sight and taking account 
of extinction due to the line of sight optical depth. This spatial intensity distribution at 
any selected wavelengths is convolved with the relevant instrumental beam profile (PSF) for 
direct comparison with observations. 
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Before applying the scheme developed above to any astrophysical source it is necessary 
to verify its reliability and quantify its accuracy. For this, we simulate the case of a single 
embedded source by "dimming" one of the two sources thereby keeping our original code in- 
tact during the test runs. This simulates a single exciting source embedded on the symmetry 
axis of the uniform density cylindrically shaped cloud. The size of the cylindrical cloud has 
been chosen such that the results from other codes using spherically symmetric geometry 
could be compared effectively. We have used the well established code CSDUST3 (Egan, 
Leung & Spagna 1988) for such a comparison. 

The radius of the cylindrical cloud used for comparing our code with the spherically 
symmetric code was identical to that of the "equivalent" spherical cloud. The length of the 
cylinder (along z-axis) is twice this radius. Figure 2 is a schematic of the cylindrical as well 
as the spherical models. 

The assumed radius of the cylinder, as well as the sphere, was 2 pc. The radius of the 
dust free cavity near the source was adopted to be .01 pc. In order to make the comparison 
possible and effective, all model parameters were made identical for both the codes. These 
include : dust size distribution; dust number density (and hence total optical depth along 
the line of sight between the embedded source and the distant observer). A typical dust 
composition has been assumed consisting of 50 % Graphite and 50 % Astronomical Silicate. 
The embedded energy source was assumed to be a single ZAMS 05 star with luminosity 
9.0 x 1O 5 L . Some parameters were explored (e.g. optical depth) to study regions of validity 
with specified accuracy, by changing them identically for both the schemes. A large range 
in the optical depth Tioo, viz., 1.0 x 10~ 3 to 0.1 was covered during the test runs. 

The comparison of the emergent spectral energy distributions (SEDs) from both the 
schemes for various optical depths are shown in Fig. 3. It may be seen that the SEDs 
match quantitatively over a wide range of optical depths from mid infrared to millimeter 
wavelengths. There are some differences at near infrared wavelengths; this is to be expected 
from the differences in the cell sizes and the geometry. However, since the main motive of 
the present study is to interpret measurements in the wavebands beyond the mid-IR, our 
code may be considered to be satisfactory. From this comparative study we conclude that 
our code is accurate upto an optical depth corresponding to tioo ~ 0.06 (ti Mm ~ 6), for the 
present choice of grid points. This limit already corresponds to much denser clouds than 
generally found in the Galactic star forming regions. 
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3 IRAS 19181+1349 



The Galactic star forming region IRAS 19181 + 1349 is a IRAS Point Source Catalog (IRAS 
PSC) source associated with the radio source G48. 60+00. The presence of a radio continuum 
suggests ongoing high mass star formation in the region. This source has been resolved in 
two components in the 210 fim map (Fig. 4b) generated from the observations using the 
TIFR 1-meter balloon borne telescope (Karnik et al, 1999). Although a single IRAS PSC 
source is associated with this star forming region, the HIRES processing of the IRAS survey 
data has led to the resolution of these two sources in the 12 & 25 /iin band maps. The map 
at 12 jum is shown in Fig. 4a. We denote FIR peak with higher R.A. as SI and that of 
lower R.A. S2, respectively, hereafter. Zoonematkermani et al. (1990) have reported four 
compact radio sources 48.603+0.026,48.609+0.027,48.606+0.023 and 48.592+0.044 in this 
region. First three of these closely match SI (positionally) while the fourth one matches 
S2. Kurtz, Churchwell & Wood (1994) have also reported three ultracompact HII regions 
which are positionally close to SI, indicating ongoing high mass star formation in a very 
early evolutionary phase. This source has also been observed using ISOCAM instrument 
onboard Infrared Space Observatory (ISO) using seven filters centred at four PAH features 
and 3 continuua (3.30, 3.72, 6.00, 6.75, 7.75, 9.63, 11.37 //m; Verma et al 1999). The 
ISOCAM images at all the bands show two extended prominent complexes with multiple 
peaks. Hence from radio as well as mid to far infrared observations the nature of IRAS 
19181+1349 is established to be of double embedded source/cluster type. 

This seems to be an ideal astrophysical example of interstellar cloud with two embedded 
sources. We attempt to model this source with our scheme to extract important physical 
parameters about the embedded energy sources as well as the intervening interstellar medium 
in IRAS 19181+1349. Next we describe the results obtained from such a study. 

The source IRAS 19181+1349 was considered as a cylindrical dust cloud with two pro- 
tostellar /ZAMS stellar sources embedded along the axis of the cylinder. The size of the 
cloud and the dust density were used as free parameters. The sum of the luminosities of 
the two embedded sources is determined by integrating the observed spectral energy dis- 
tribution (SED). The observations used for this integration includes the four IRAS bands 
and the two TIFR bands. This total luminosity is treated as an observational constraint in 
the modelling. Further observational constraints include : the shape of the SED which was 
obtained from HIRES-IRAS, TIFR and ISO observations, and the structural morphology of 
IRAS 19181+1349 as reflected by the isophote contours of the high angular resolution maps 
at 12 & 210 /iin. As our code does not deal with the physics of PAH emission at present, 
only the data in ISOCAM filters sampling the continuua have been used to constrain the 
model. 
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Two completely independent approaches have been followed in modelling IRAS 19181+1349 
using our scheme. In each approach, all the model parameters are floated to obtain the best 
fit model. The parameters fine tuned to achieve the best fit model are : luminosities of indi- 
vidual embedded sources; geometrical size details of the cylindrical cloud (including the size 
of the cavity); and the dust density / optical depth. The main aim of the first approach is to 
optimize the fit to the observed SED (hereafter Msed)- The second approach optimizes the 
one dimensional radial intensity profiles at selected wavelengths covering mid to far infrared 
bands (hereafter M RC ). The radial profiles are taken along geometrically interesting axes, 
viz., along the line joining the two embedded sources and along lines perpendicular to the 
cylinder axis and passing through either of the two sources (see Fig. 4). Whereas the first 
approach gives precedence to the overall energetics the latter gives more importance to the 
structural details in the isophotes, particularly close to the embedded sources and the region 
between them. The actual reality may lie somewhere in between these two approaches. 

The axes where radial cuts have been taken are also displayed on the 12 & 210 fxm maps 
in Fig. 4. The fits to the observed SED for both Msed as well as M R c are presented in Fig. 
5. Comparison of the radial cuts at 12 & 210 /mi along the three axes between the observed 
maps and the best fit M RC model are shown in Fig. 6. 

The Msed approach fits the observed SED very well right through near-IR to sub- mm, 
though there is some discrepancy in the far-IR (Fig. 5). On the other hand, fit from the 
Mrc approach is resonable for A > 25 /zm only. The radial cuts at 210 /im fit the observed 
data very well. However, at 12 /im although the model predictions qualitatively agree with 
the data, there are discrepancies viz., the extent of emission along vertical axis and the 
relative contrast of the minima between the two sources. From the above it appears that 
the M RC approach is very sensitive to the exact distribution of sources, specially at shorter 
wavelengths (which traces much hotter dust and hence physically closer to the exciting 
source). One possibility for these discrepancies is that our model is too simplified than the 
actual reality in IRAS 19181+1349. For example there may be a distribution of sources 
along the line joining SI and S2 and away from it, but with more concentration near SI and 
S2. This way the discrepancy for M R c in the mid-IR part of SED will also be explained. 
Assuming a constant dust density in the cloud may also be a major reason for the discrepancy. 

The results of modelling IRAS 19181+1349 in the form of best fit model parameters as 
found under both the approaches {Msed and M R c) are presented in Table 1. The value of 
t 10 o found from both the approaches are well within the range of validity of the code, i.e. 
much less than 0.06. The luminosities of both the sources are similar and the dust is found to 
be predominantly (80%) Silicate with no Silicon Carbide, the rest (20%) being Graphite. It 
is remarkable that almost all the important parameters are identical as found from both the 



9 



approaches. The only difference is seen at the physical sizes of the cavities. The remarkable 
similar physical parameters obtained from both the approaches gives good confidence on the 
derived astrophysical parameters for the Galactic star forming region IRAS 19181 + 1349. 

This scheme will be useful to model other similar double sources. A huge database from 
ISO (SEDs as well as images leading to radial cuts, covering 2.5 to 200 fxm of the spectrum) 
will become available soon for such modelling. 

4 Summary 

High angular resolution mid / far infrared maps of many Galactic star forming regions show 
evidence for multiple embedded energy sources in the corresponding interstellar clouds. With 
the aim of studying such star forming regions with two embedded sources, a scheme has 
been developed to carry out radiative transfer calculations in a uniform density dust cloud 
of cylindrical geometry, and which includes the effect of isotropic scattering in addition to 
the absorption and emission processes. In addition to the luminosities of the two embedded 
energy sources, the cylindrical cloud size, separation between the two sources, dust density 
and composition are the parameters for the modelling. The accuracy of our scheme has been 
tested by comparing the results with a well established 1-D code. 

An attempt was made to model the Galactic star forming region associated with IRAS 
19181+1349 which shows double peaks in the mid and far infrared maps using the scheme 
described here. Two independent approaches were employed to find the best fit models for 
IRAS 19181+1349. Whereas the first (Msed) approach aims to fit the observed SED best; 
the latter (Mac) aims to fit the radial intensity profiles (along a few important axes) at mid to 
far infrared wavebands. Interestingly most of the crucial model parameters like luminosities, 
effective temperatures, dust composition, optical depth etc turn out to be identical under 
both the approaches. 
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Table 1: Best fit model Parameters for M RC and M S ed schemes. 



Parameter M RC 



^SED 



R c (pc)(Sourcel) 0.11 
R c (pc)(Source2) 0.09 



0.03 
0.01 



Lsourcel (L© ) 6.2 XlO 5 

L S ource2 (L© ) 6.3 XlO 5 



6.2 xlO 5 

6.3 xlO 5 



TlOOcyl (/PC) 0.024 
Rcyl (PC) 2.1 



0.028 
1.8 



Graphite:Silicate:SiC 20:80:0 



20:80:0 
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Figure 1: The geometry of the cylindrical cloud considered in the present study. The 
embedded energy sources lie along the axis of the cylinder (at the centers of the two small 
cylinders) seperated by distance D. The two small cylinders represent the dust free cavities 
near the respective sources. The z-axis is defined to be along the symmetry axis of the 
cylinder. The cloud is divided into n z disks along z-axis, and each disk into n r concentric 
rings. 
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Figure 3: The comparison of predicted SEDs from the two schemes. The continuous lines 
represent CSDUST3 while the dashed lines refer to our scheme. The SEDs correspond to 
r 100 of - (a) 0.001, (b) .01, (c) .06 and (d) 0.1. 
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Figure 4: The infrared maps of IRAS 19181+1349 for (a) 12 /j,m (b) 210 /mi. Contours are 
drawn at 10, 20, 30, 40, 50, 60, 70, 80, 90, 95 percent of peak intensities 26 and 303 Jy/ (arc 
min) 2 respectively. The straight lines show the positions of radial and axial cuts compared 
with the model in Fig. || 



15 



10000 




log (A fim ) 



Figure 5: The SEDs for the source IRAS 19181 + 1349 predicted by the radiative trans- 
fer models. The solid and the dashed lines represent best fit models from the M$ed and 
Mrc approaches, respectively. The crosses represent HIRES (processed IRAS), TIFR and 
ISOCAM continuum flux densities. The solid triangles represent ISOCAM flux densities 
measured through filters centered on the PAH features. 
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Figure 6: The axial and radial cuts predicted by the radiative transfer model (M RC ) 
(continuous lines), compared with the observations (crosses), at 12 and 210 /im. The spatial 
positions of the cuts are displayed in Fig. £|. 
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